/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains classes for grain size distributions.

// $Id$

#ifndef __grain_size__
#define __grain_size__

#include "blitz/array.h"
#include "mcrx-types.h"
#include "mcrx-units.h"
#include "constants.h"

namespace mcrx {
  class grain_size_distribution;
  class wd01_distribution;
  class wd01_graphite_distribution;
  class wd01_silicate_distribution;
};


/** Abstract base class for grain size distributions.  Only declares
    the virtual function dn_da.  */ 
class mcrx::grain_size_distribution {
public:
  /** This function returns dn/(da*M_dust), that is the number of dust
      grains per size interval, normalized to 1 M_sun of dust. */
  virtual T_float dn_da(T_float a) const=0;
};


/** Class representing size distributions of the functional form used
    by Weingartner and Draine 2001. The derived classes define the
    parameters for the distribution. The units of mass and size are
    taken from the unit map supplied to the constructor. */ 
class mcrx::wd01_distribution : public mcrx::grain_size_distribution {
public:
  T_float alpha_;
  T_float beta_;
  T_float a_t_;      ///< Has units of length.
  T_float a_c_;      ///< Has units of length.
  T_float C_;
  T_float massnorm_; ///< Mass normalization. Number of H atoms per unit mass.
  T_float sizenorm_; ///< Size normalization. Length of size unit in m.

  // graphite-only numbers -- only defined for graphite grains
  T_float m_C_;      ///< Carbon atom mass in kg (1.99e-26).
  T_float rho_;      ///< Carbon density in kg/m3 (2.24e3).
  T_float b_C1_;
  T_float b_C2_;
  T_float a0_1_;
  T_float a0_2_;
  T_float sigma_1_;
  T_float sigma_2_;

  T_float F(T_float a, T_float beta, T_float a_t) const {
    return (beta>=0) ? 1 + beta*a/a_t : 1./(1 - beta*a/a_t);};

  T_float B(T_float a0, T_float b_C, T_float sigma) const {
    return 3/pow(2*constants::pi, 1.5)*
      exp(-4.5*sigma*sigma)/(rho_*a0*a0*a0*sigma)*
      b_C*m_C_/(1 + erf(3*sigma/sqrt(2.)+log(a0/3.5e-10)/(sigma*sqrt(2.))));};

  T_float D(T_float a) const {
    return 
      B(a0_1_,b_C1_, sigma_1_)/a*exp(-0.5*pow(log(a/a0_1_)/sigma_1_, 2)) +
      B(a0_2_,b_C2_, sigma_2_)/a*exp(-0.5*pow(log(a/a0_2_)/sigma_2_, 2));
  };

public:
  wd01_distribution(const T_unit_map& units) : 
   alpha_(0), beta_(0), a_t_(0), a_c_(0), C_(0), massnorm_(0),
   sizenorm_(0), m_C_(1.99447e-26), rho_(2240.), b_C1_(0), b_C2_(0), 
   a0_1_(1), a0_2_(1), sigma_1_(1), sigma_2_(1) 
  {
    // 1.0973e59 is the number of H atoms per Msun of dust in the WD01
    // model, so we need to convert the normalization to whatever
    // units we are using
    massnorm_ = 1.0973e59* units::convert(units.get("mass"), "Msun");
    sizenorm_ = units::convert(units.get("length"), "m");
  };

  /** This function returns dn/(da*M_dust), that is the number of dust
      grains per size interval, normalized to unit dust mass, for the
      WD01 dust model. */
  virtual T_float dn_da(T_float a) const {
    a *= sizenorm_;
    // because we return dn/da, we need to convert the return value to
    // reflect the unit of da in the denominator, too.
    const T_float dnda = massnorm_ * sizenorm_ *
      (D(a) + C_/a*pow(a/a_t_,alpha_)*F(a, beta_, a_t_)*
       ((a<a_t_)? 1 : exp(-pow((a-a_t_)/a_c_,3))));
    assert(dnda==dnda);
    return dnda;
  };
  
  BZ_DECLARE_MEMBER_FUNCTION(wd01_distribution, dn_da);
};


/** This class is a wd01_distribution with the particular parameters
    describing the Weingartner and Draine graphite distribution. The
    model string specifies which of the particular models are used. */ 
class mcrx::wd01_graphite_distribution : public mcrx::wd01_distribution {
public:
  typedef blitz::TinyVector<T_float, 10> params;

private:
  /** Static map of parameter set name to the parameter vector that
      goes into the size distribution. Used by the constructor to set
      the size distribution parameters. */
  static std::map<std::string, params> paramsets;

  /** Set the size distribution parameters according to the paramset
      object */
  void set_params(const params& p);

public:
  /** Create a WD01 graphite size distribution with the parameters of
      the secified model (from WD01 paper). */
  wd01_graphite_distribution (const std::string& model,
			      const T_unit_map& units);
};


/** This class is a wd01_distribution with the particular parameters
    describing the Weingartner and Draine silicate distribution. */ 
class mcrx::wd01_silicate_distribution : public mcrx::wd01_distribution {
public:
  typedef blitz::TinyVector<T_float, 5> params;

private:
  /** Static map of parameter set name to the parameter vector that
      goes into the size distribution. Used by the constructor to set
      the size distribution parameters. */
  static std::map<std::string, params> paramsets;

  /** Set the size distribution parameters according to the paramset
      object */
  void set_params(const params& p);

public:
  /** Create a WD01 silicate size distribution with the parameters of
      the secified model (from WD01 paper). */
  wd01_silicate_distribution (const std::string& model,
			      const T_unit_map& units);
};

#endif
